Model Spesifications

# Essential libraries
library(corrplot)
library(ggeffects)
library(ggplot2)
library(DHARMa)
library(glmmTMB)
library(dplyr)
library(lubridate)
library(emmeans)
library(tidyverse)
library(performance)

1. Data Import and Preparation

df <- read.csv("./data/Mock_Mosquito_MixedEffects_w_slopes_Dataset.csv")
df$Trap_type = factor(df$Trap_type)
df$Mosquito_species = factor(df$Mosquito_species)

head(df)

2. Model Specification

Defining Random Effects

We start by first defining the random effects part. There are two random effects defined in our model, one for the nested structure in our data and the other for the mosquito species. By defining random intercepts in m1 model, we assume that different cities have different baseline abundance levels, different sites within the same city have different baseline abundance, and different mosquito species have different baseline abundance levels. However, the slopes remain consistent meaning all cities and species respond to green coverage the same way.

m1_baseline <- glmmTMB(
  Mosquito_abundance ~ Green_area_coverage * Trap_type +
                       Temperature * Relative_humidity +
                       Host_abundance * Predator_abundance +
                       Blue_area_coverage +
                       factor(Month) +
                       (1 | City / Location.Site) +      # nested random effects
                       (1 | Mosquito_species),           # random intercept by species
  family = nbinom2(),
  data = df,
  REML = TRUE
)

summary(m1_baseline)

 Family: nbinom2  ( log )
Formula:          Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature *  
    Relative_humidity + Host_abundance * Predator_abundance +  
    Blue_area_coverage + factor(Month) + (1 | City/Location.Site) +      (1 | Mosquito_species)
Data: df

      AIC       BIC    logLik -2*log(L)  df.resid 
  18366.2   18515.3   -9157.1   18314.2      2262 

Random effects:

Conditional model:
 Groups             Name        Variance Std.Dev.
 Location.Site:City (Intercept) 0.12116  0.3481  
 City               (Intercept) 0.12336  0.3512  
 Mosquito_species   (Intercept) 0.05108  0.2260  
Number of obs: 2288, groups:  Location.Site:City, 54; City, 6; Mosquito_species, 5

Dispersion parameter for nbinom2 family (): 1.98 

Conditional model:
                                         Estimate Std. Error z value Pr(>|z|)    
(Intercept)                            -1.867e+00  7.967e-01  -2.344 0.019084 *  
Green_area_coverage                     8.413e-03  2.284e-03   3.683 0.000231 ***
Trap_typeLight_CO2                      1.175e-01  7.294e-02   1.611 0.107204    
Temperature                             7.078e-02  3.399e-02   2.082 0.037300 *  
Relative_humidity                       3.435e-02  8.128e-03   4.226 2.38e-05 ***
Host_abundance                          2.017e-02  9.056e-03   2.227 0.025963 *  
Predator_abundance                      2.250e-03  6.004e-03   0.375 0.707912    
Blue_area_coverage                      6.065e-03  3.833e-03   1.583 0.113513    
factor(Month)2                         -5.252e-02  8.426e-02  -0.623 0.533063    
factor(Month)3                          1.012e-01  9.255e-02   1.093 0.274366    
factor(Month)4                          1.588e-01  1.064e-01   1.492 0.135757    
factor(Month)5                          3.574e-01  1.201e-01   2.975 0.002934 ** 
factor(Month)6                          2.764e-01  1.308e-01   2.113 0.034565 *  
factor(Month)7                          3.973e-01  1.215e-01   3.270 0.001075 ** 
factor(Month)8                          1.085e-01  1.059e-01   1.025 0.305318    
factor(Month)9                          1.525e-01  9.086e-02   1.678 0.093279 .  
factor(Month)10                         8.593e-02  8.213e-02   1.046 0.295435    
factor(Month)11                         3.149e-04  7.928e-02   0.004 0.996830    
factor(Month)12                        -1.800e-03  8.171e-02  -0.022 0.982421    
Green_area_coverage:Trap_typeLight_CO2  5.819e-03  1.265e-03   4.599 4.26e-06 ***
Temperature:Relative_humidity          -8.019e-04  4.515e-04  -1.776 0.075746 .  
Host_abundance:Predator_abundance       1.507e-05  1.160e-04   0.130 0.896639    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

baseline_aic <- AIC(m1_baseline)
cat("Baseline AIC:", baseline_aic, "\n")

Baseline AIC: 18366.19 

In the second model, we consider random slopes for green levels in the cities, meaning that the effect of urban green differs per city.

m2_city_green <- glmmTMB(
  Mosquito_abundance ~ Green_area_coverage * Trap_type +
                       Temperature * Relative_humidity +
                       Host_abundance * Predator_abundance +
                       Blue_area_coverage +
                       factor(Month) +
                      (1 + Green_area_coverage | City) +          # City random intercept + slope
                      (1 | Location.Site:City) +        
                      (1 | Mosquito_species),
  family = nbinom2(),
  data = df,
  REML =TRUE
)


# Check convergence
check_convergence(m2_city_green)

[1] TRUE

m2_city_green$sdr$pdHess

[1] TRUE

cat("AIC difference:", AIC(m2_city_green) - baseline_aic, "\n")

AIC difference: -2.716824 

lrt_m2 <- anova(m1_baseline, m2_city_green)$Chisq[2]

p_m2 <- 0.5 * pchisq(lrt_m2, df = 2, lower.tail = FALSE) + 
           0.5 * pchisq(lrt_m2, df = 1, lower.tail = FALSE)
cat("p-value:", p_m2, "\n")

p-value: 0.02217061 

AIC difference is negative and the \(H_0: \text{Cities do not differ in how urban green level affects mosquito abundance}\) is rejected (p<0.05), the random slope of urban green level is significant.

m3_city_temp <- glmmTMB(
  Mosquito_abundance ~ Green_area_coverage * Trap_type +
                       Temperature * Relative_humidity +
                       Host_abundance * Predator_abundance +
                       Blue_area_coverage +
                       factor(Month) +
                      (1 + Temperature | City) +      
                      (1 | Location.Site:City) +        
                      (1 | Mosquito_species),
  family = nbinom2(),
  data = df,
  REML = TRUE
)

# Check convergence
check_convergence(m3_city_temp)

[1] FALSE

m3_city_temp$sdr$pdHess

[1] TRUE

cat("AIC difference:", AIC(m3_city_temp) - baseline_aic, "\n")

AIC difference: 3.928898 

lrt_m3 <- anova(m1_baseline, m3_city_temp)$Chisq[2]

p_m3 <- 0.5 * pchisq(lrt_m3, df = 2, lower.tail = FALSE) +
        0.5 * pchisq(lrt_m3, df = 1, lower.tail = FALSE)
cat("p-value for Temp slope:", p_m3, "\n")

p-value for Temp slope: 0.8774065 

AIC difference is positive and the \(H_0: \text{Cities do not differ in how temperature affects mosquito abundance}\) is retained (p>0.05), the random slope of temperature is not significant.

m4_both_slopes <- glmmTMB(
  Mosquito_abundance ~ Green_area_coverage * Trap_type +
                       Temperature * Relative_humidity +
                       Host_abundance * Predator_abundance +
                       Blue_area_coverage +
                       factor(Month) +
                       (1 + Green_area_coverage + Temperature | City) +
                       (1 | Location.Site:City) +
                       (1 | Mosquito_species),
  family = nbinom2(),
  data = df,
  REML = TRUE
)

# Check convergence
check_convergence(m4_both_slopes)

[1] FALSE

m4_both_slopes$sdr$pdHess

[1] TRUE

# Compare m2 -> m4 (does adding temperature slope to m2 improve fit?)
lrt_m2_m4 <- anova(m2_city_green, m4_both_slopes)$Chisq[2]
cat("LRT stat (m2 -> m4):", lrt_m2_m4, "\n")

LRT stat (m2 -> m4): 0.204401 

p_m2_m4 <- 0.5 * pchisq(lrt_m2_m4, df = 2, lower.tail = FALSE) + 
           0.5 * pchisq(lrt_m2_m4, df = 1, lower.tail = FALSE)
cat("p-value (m2 -> m4):", p_m2_m4, "\n")

p-value (m2 -> m4): 0.7770201 

# Compare m3 -> m4 (does adding Green slope to m3 improve fit?)
lrt_m3_m4 <- anova(m3_city_temp, m4_both_slopes)$Chisq[2]
cat("LRT stat (m3 -> m4):", lrt_m3_m4, "\n")

LRT stat (m3 -> m4): 6.850123 

p_m3_m4 <- 0.5 * pchisq(lrt_m3_m4, df = 2, lower.tail = FALSE) + 
           0.5 * pchisq(lrt_m3_m4, df = 1, lower.tail = FALSE)
cat("p-value (m3 -> m4):", p_m3_m4, "\n")

p-value (m3 -> m4): 0.0207054 

# Compare all models by AIC 
models <- list(m1_baseline = m1_baseline, m2 = m2_city_green, m3 = m3_city_temp, m4 = m4_both_slopes)
aic_tab <- data.frame(
  model = names(models),
  AIC = sapply(models, AIC),
  pdHess = sapply(models, function(m) { if (is.null(m$sdr)) NA else m$sdr$pdHess })
)
aic_tab <- aic_tab %>% arrange(AIC)
print(aic_tab)

Across the models tested, the lowest AIC value was tested was achieved by model 2, for the remainder of the model selection procedure we will use the random effects defined in model 2. Now, we will focus on getting a more parsimonious fixed effects part by removing interaction terms in a step wise manner.

m_full <- glmmTMB(
  Mosquito_abundance ~ Green_area_coverage * Trap_type +
                       Temperature * Relative_humidity +
                       Host_abundance * Predator_abundance +
                       Blue_area_coverage +
                       factor(Month) +
                       (1 + Green_area_coverage | City) +          # City random intercept + slope
                      (1 | Location.Site:City) +        
                      (1 | Mosquito_species),
  family = nbinom2(),
  data = df,
  REML = FALSE  # ML for fixed effects comparison
)

summary(m_full)

 Family: nbinom2  ( log )
Formula:          Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature *  
    Relative_humidity + Host_abundance * Predator_abundance +  
    Blue_area_coverage + factor(Month) + (1 + Green_area_coverage |  
    City) + (1 | Location.Site:City) + (1 | Mosquito_species)
Data: df

      AIC       BIC    logLik -2*log(L)  df.resid 
  18216.1   18376.7   -9080.0   18160.1      2260 

Random effects:

Conditional model:
 Groups             Name                Variance  Std.Dev. Corr  
 City               (Intercept)         6.083e-02 0.246641       
                    Green_area_coverage 3.707e-05 0.006088 -0.37 
 Location.Site:City (Intercept)         8.948e-02 0.299125       
 Mosquito_species   (Intercept)         4.596e-02 0.214383       
Number of obs: 2288, groups:  City, 6; Location.Site:City, 54; Mosquito_species, 5

Dispersion parameter for nbinom2 family ():    2 

Conditional model:
                                         Estimate Std. Error z value Pr(>|z|)    
(Intercept)                            -1.862e+00  7.828e-01  -2.379 0.017359 *  
Green_area_coverage                     7.260e-03  3.286e-03   2.209 0.027153 *  
Trap_typeLight_CO2                      1.245e-01  7.270e-02   1.713 0.086802 .  
Temperature                             6.912e-02  3.385e-02   2.042 0.041182 *  
Relative_humidity                       3.430e-02  8.095e-03   4.237 2.27e-05 ***
Host_abundance                          2.047e-02  8.994e-03   2.276 0.022868 *  
Predator_abundance                      2.477e-03  5.961e-03   0.415 0.677825    
Blue_area_coverage                      7.381e-03  3.584e-03   2.060 0.039431 *  
factor(Month)2                         -4.997e-02  8.398e-02  -0.595 0.551836    
factor(Month)3                          1.074e-01  9.227e-02   1.164 0.244285    
factor(Month)4                          1.671e-01  1.062e-01   1.574 0.115443    
factor(Month)5                          3.700e-01  1.198e-01   3.088 0.002016 ** 
factor(Month)6                          2.886e-01  1.304e-01   2.213 0.026866 *  
factor(Month)7                          4.082e-01  1.212e-01   3.368 0.000757 ***
factor(Month)8                          1.170e-01  1.056e-01   1.108 0.268040    
factor(Month)9                          1.577e-01  9.056e-02   1.741 0.081608 .  
factor(Month)10                         8.754e-02  8.185e-02   1.070 0.284832    
factor(Month)11                         5.986e-04  7.900e-02   0.008 0.993954    
factor(Month)12                        -1.046e-03  8.141e-02  -0.013 0.989752    
Green_area_coverage:Trap_typeLight_CO2  5.691e-03  1.262e-03   4.511 6.45e-06 ***
Temperature:Relative_humidity          -7.912e-04  4.496e-04  -1.760 0.078454 .  
Host_abundance:Predator_abundance       1.126e-05  1.152e-04   0.098 0.922116    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

We start by removing the least significant interaction term which is Host_abundance:Predator_abundance .

m_wo_host_pre <- glmmTMB(
  Mosquito_abundance ~ Green_area_coverage * Trap_type +
                       Temperature * Relative_humidity +
                       Host_abundance + Predator_abundance +
                       Blue_area_coverage +
                       factor(Month) +
                       (1 + Green_area_coverage | City) +          # City random intercept + slope
                      (1 | Location.Site:City) +        
                      (1 | Mosquito_species),
  family = nbinom2(),
  data = df,
  REML = FALSE  # ML for fixed effects comparison
)

summary(m_wo_host_pre)

 Family: nbinom2  ( log )
Formula:          Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature *  
    Relative_humidity + Host_abundance + Predator_abundance +  
    Blue_area_coverage + factor(Month) + (1 + Green_area_coverage |  
    City) + (1 | Location.Site:City) + (1 | Mosquito_species)
Data: df

      AIC       BIC    logLik -2*log(L)  df.resid 
  18214.1   18369.0   -9080.0   18160.1      2261 

Random effects:

Conditional model:
 Groups             Name                Variance  Std.Dev. Corr  
 City               (Intercept)         6.066e-02 0.246290       
                    Green_area_coverage 3.711e-05 0.006092 -0.37 
 Location.Site:City (Intercept)         8.961e-02 0.299345       
 Mosquito_species   (Intercept)         4.597e-02 0.214411       
Number of obs: 2288, groups:  City, 6; Location.Site:City, 54; Mosquito_species, 5

Dispersion parameter for nbinom2 family ():    2 

Conditional model:
                                         Estimate Std. Error z value Pr(>|z|)    
(Intercept)                            -1.9033561  0.6709093  -2.837 0.004554 ** 
Green_area_coverage                     0.0072419  0.0032831   2.206 0.027399 *  
Trap_typeLight_CO2                      0.1243781  0.0726983   1.711 0.087103 .  
Temperature                             0.0690642  0.0349110   1.978 0.047896 *  
Relative_humidity                       0.0342986  0.0082954   4.135 3.56e-05 ***
Host_abundance                          0.0213223  0.0021267  10.026  < 2e-16 ***
Predator_abundance                      0.0030338  0.0017445   1.739 0.082024 .  
Blue_area_coverage                      0.0074033  0.0035795   2.068 0.038618 *  
factor(Month)2                         -0.0501047  0.0839359  -0.597 0.550547    
factor(Month)3                          0.1074290  0.0921659   1.166 0.243775    
factor(Month)4                          0.1675267  0.1060089   1.580 0.114037    
factor(Month)5                          0.3701742  0.1199362   3.086 0.002026 ** 
factor(Month)6                          0.2887538  0.1304774   2.213 0.026894 *  
factor(Month)7                          0.4082087  0.1213271   3.365 0.000767 ***
factor(Month)8                          0.1170888  0.1055235   1.110 0.267172    
factor(Month)9                          0.1578000  0.0904607   1.744 0.081089 .  
factor(Month)10                         0.0876332  0.0818319   1.071 0.284218    
factor(Month)11                         0.0006265  0.0790036   0.008 0.993673    
factor(Month)12                        -0.0008575  0.0813868  -0.011 0.991593    
Green_area_coverage:Trap_typeLight_CO2  0.0056932  0.0012615   4.513 6.38e-06 ***
Temperature:Relative_humidity          -0.0007906  0.0004658  -1.697 0.089621 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

anova(m_full, m_wo_host_pre)

Data: df
Models:
m_wo_host_pre: Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature * , zi=~0, disp=~1
m_wo_host_pre:     Relative_humidity + Host_abundance + Predator_abundance + , zi=~0, disp=~1
m_wo_host_pre:     Blue_area_coverage + factor(Month) + (1 + Green_area_coverage | , zi=~0, disp=~1
m_wo_host_pre:     City) + (1 | Location.Site:City) + (1 | Mosquito_species), zi=~0, disp=~1
m_full: Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature * , zi=~0, disp=~1
m_full:     Relative_humidity + Host_abundance * Predator_abundance + , zi=~0, disp=~1
m_full:     Blue_area_coverage + factor(Month) + (1 + Green_area_coverage | , zi=~0, disp=~1
m_full:     City) + (1 | Location.Site:City) + (1 | Mosquito_species), zi=~0, disp=~1
              Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
m_wo_host_pre 27 18214 18369  -9080    18160                         
m_full        28 18216 18377  -9080    18160 0.0095      1     0.9222

Removing the Host_abundance:Predator_abundance interaction does not significantly worsen model fit, then the interaction is not needed. We move on to the next interaction.

m_wo_temp_hum <- glmmTMB(
  Mosquito_abundance ~ Green_area_coverage * Trap_type +
                       Temperature + Relative_humidity +
                       Host_abundance + Predator_abundance +
                       Blue_area_coverage +
                       factor(Month) +
                       (1 + Green_area_coverage | City) +          # City random intercept + slope
                      (1 | Location.Site:City) +        
                      (1 | Mosquito_species),
  family = nbinom2(),
  data = df,
  REML = FALSE  # ML for fixed effects comparison
)

summary(m_wo_temp_hum)

 Family: nbinom2  ( log )
Formula:          Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature +  
    Relative_humidity + Host_abundance + Predator_abundance +  
    Blue_area_coverage + factor(Month) + (1 + Green_area_coverage |  
    City) + (1 | Location.Site:City) + (1 | Mosquito_species)
Data: df

      AIC       BIC    logLik -2*log(L)  df.resid 
  18215.2   18364.3   -9081.6   18163.2      2262 

Random effects:

Conditional model:
 Groups             Name                Variance  Std.Dev. Corr  
 City               (Intercept)         6.149e-02 0.247970       
                    Green_area_coverage 3.759e-05 0.006131 -0.38 
 Location.Site:City (Intercept)         8.909e-02 0.298473       
 Mosquito_species   (Intercept)         4.585e-02 0.214137       
Number of obs: 2288, groups:  City, 6; Location.Site:City, 54; Mosquito_species, 5

Dispersion parameter for nbinom2 family ():    2 

Conditional model:
                                         Estimate Std. Error z value Pr(>|z|)    
(Intercept)                            -0.9280810  0.3606140  -2.574 0.010064 *  
Green_area_coverage                     0.0071653  0.0032923   2.176 0.029527 *  
Trap_typeLight_CO2                      0.1242485  0.0727127   1.709 0.087496 .  
Temperature                             0.0113019  0.0082549   1.369 0.170962    
Relative_humidity                       0.0211927  0.0031634   6.699 2.09e-11 ***
Host_abundance                          0.0212425  0.0021284   9.981  < 2e-16 ***
Predator_abundance                      0.0031070  0.0017451   1.780 0.075007 .  
Blue_area_coverage                      0.0074403  0.0035686   2.085 0.037073 *  
factor(Month)2                         -0.0599509  0.0838062  -0.715 0.474392    
factor(Month)3                          0.0958339  0.0921129   1.040 0.298156    
factor(Month)4                          0.1601716  0.1060906   1.510 0.131104    
factor(Month)5                          0.3753467  0.1198993   3.131 0.001745 ** 
factor(Month)6                          0.2967466  0.1303588   2.276 0.022823 *  
factor(Month)7                          0.4141824  0.1212072   3.417 0.000633 ***
factor(Month)8                          0.1090555  0.1055911   1.033 0.301693    
factor(Month)9                          0.1453758  0.0903559   1.609 0.107633    
factor(Month)10                         0.0820182  0.0818185   1.002 0.316131    
factor(Month)11                         0.0004793  0.0790249   0.006 0.995160    
factor(Month)12                        -0.0005138  0.0814086  -0.006 0.994964    
Green_area_coverage:Trap_typeLight_CO2  0.0057447  0.0012610   4.556 5.22e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

anova(m_wo_host_pre, m_wo_temp_hum)

Data: df
Models:
m_wo_temp_hum: Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature + , zi=~0, disp=~1
m_wo_temp_hum:     Relative_humidity + Host_abundance + Predator_abundance + , zi=~0, disp=~1
m_wo_temp_hum:     Blue_area_coverage + factor(Month) + (1 + Green_area_coverage | , zi=~0, disp=~1
m_wo_temp_hum:     City) + (1 | Location.Site:City) + (1 | Mosquito_species), zi=~0, disp=~1
m_wo_host_pre: Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature * , zi=~0, disp=~1
m_wo_host_pre:     Relative_humidity + Host_abundance + Predator_abundance + , zi=~0, disp=~1
m_wo_host_pre:     Blue_area_coverage + factor(Month) + (1 + Green_area_coverage | , zi=~0, disp=~1
m_wo_host_pre:     City) + (1 | Location.Site:City) + (1 | Mosquito_species), zi=~0, disp=~1
              Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
m_wo_temp_hum 26 18215 18364 -9081.6    18163                           
m_wo_host_pre 27 18214 18369 -9080.0    18160 3.1134      1    0.07765 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Removing the Temp x Humidity interaction does not worsen model fit, the interaction should be removed. The remaining interaction is significant, therefore we move to the main effects.

m_wo_blue <- glmmTMB(
  Mosquito_abundance ~ Green_area_coverage * Trap_type +
                       Temperature + Relative_humidity +
                       Host_abundance + Predator_abundance +
                       factor(Month) +  
                       (1 + Green_area_coverage | City) +          # City random intercept + slope
                      (1 | Location.Site:City) +        
                      (1 | Mosquito_species),
  family = nbinom2(),
  data = df,
  REML = FALSE
)

summary(m_wo_blue)

 Family: nbinom2  ( log )
Formula:          Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature +  
    Relative_humidity + Host_abundance + Predator_abundance +  
    factor(Month) + (1 + Green_area_coverage | City) + (1 | Location.Site:City) +      (1 | Mosquito_species)
Data: df

      AIC       BIC    logLik -2*log(L)  df.resid 
  18217.4   18360.8   -9083.7   18167.4      2263 

Random effects:

Conditional model:
 Groups             Name                Variance  Std.Dev. Corr  
 City               (Intercept)         9.077e-02 0.30127        
                    Green_area_coverage 3.587e-05 0.00599  -0.43 
 Location.Site:City (Intercept)         9.713e-02 0.31166        
 Mosquito_species   (Intercept)         4.650e-02 0.21565        
Number of obs: 2288, groups:  City, 6; Location.Site:City, 54; Mosquito_species, 5

Dispersion parameter for nbinom2 family ():    2 

Conditional model:
                                         Estimate Std. Error z value Pr(>|z|)    
(Intercept)                            -0.7160550  0.3541022  -2.022  0.04316 *  
Green_area_coverage                     0.0074203  0.0032834   2.260  0.02382 *  
Trap_typeLight_CO2                      0.1237693  0.0726909   1.703  0.08863 .  
Temperature                             0.0112996  0.0082553   1.369  0.17107    
Relative_humidity                       0.0212518  0.0031635   6.718 1.85e-11 ***
Host_abundance                          0.0211876  0.0021288   9.953  < 2e-16 ***
Predator_abundance                      0.0030344  0.0017452   1.739  0.08208 .  
factor(Month)2                         -0.0604923  0.0838047  -0.722  0.47040    
factor(Month)3                          0.0956423  0.0921279   1.038  0.29920    
factor(Month)4                          0.1595994  0.1060774   1.505  0.13244    
factor(Month)5                          0.3750395  0.1198909   3.128  0.00176 ** 
factor(Month)6                          0.2964331  0.1303665   2.274  0.02298 *  
factor(Month)7                          0.4142905  0.1211931   3.418  0.00063 ***
factor(Month)8                          0.1085706  0.1055968   1.028  0.30387    
factor(Month)9                          0.1443938  0.0903574   1.598  0.11004    
factor(Month)10                         0.0821350  0.0818163   1.004  0.31543    
factor(Month)11                         0.0006321  0.0790305   0.008  0.99362    
factor(Month)12                        -0.0016955  0.0814145  -0.021  0.98339    
Green_area_coverage:Trap_typeLight_CO2  0.0057879  0.0012602   4.593 4.37e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

anova(m_wo_temp_hum, m_wo_blue)

Data: df
Models:
m_wo_blue: Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature + , zi=~0, disp=~1
m_wo_blue:     Relative_humidity + Host_abundance + Predator_abundance + , zi=~0, disp=~1
m_wo_blue:     factor(Month) + (1 + Green_area_coverage | City) + (1 | Location.Site:City) + , zi=~0, disp=~1
m_wo_blue:     (1 | Mosquito_species), zi=~0, disp=~1
m_wo_temp_hum: Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature + , zi=~0, disp=~1
m_wo_temp_hum:     Relative_humidity + Host_abundance + Predator_abundance + , zi=~0, disp=~1
m_wo_temp_hum:     Blue_area_coverage + factor(Month) + (1 + Green_area_coverage | , zi=~0, disp=~1
m_wo_temp_hum:     City) + (1 | Location.Site:City) + (1 | Mosquito_species), zi=~0, disp=~1
              Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
m_wo_blue     25 18217 18361 -9083.7    18167                           
m_wo_temp_hum 26 18215 18364 -9081.6    18163 4.2241      1    0.03985 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Removal of blue zone coverage worsened the model, we retain it. Month shows mixed significance, so we will try removing it and compare the model fits.

m_wo_month <- glmmTMB(
  Mosquito_abundance ~ Green_area_coverage * Trap_type +
                       Temperature + Relative_humidity +
                       Host_abundance + Predator_abundance +
                       Blue_area_coverage +
                       (1 + Green_area_coverage | City) +          # City random intercept + slope
                      (1 | Location.Site:City) +        
                      (1 | Mosquito_species),
  family = nbinom2(),
  data = df,
  REML = FALSE
)

anova(m_wo_temp_hum, m_wo_month)

Data: df
Models:
m_wo_month: Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature + , zi=~0, disp=~1
m_wo_month:     Relative_humidity + Host_abundance + Predator_abundance + , zi=~0, disp=~1
m_wo_month:     Blue_area_coverage + (1 + Green_area_coverage | City) + (1 | , zi=~0, disp=~1
m_wo_month:     Location.Site:City) + (1 | Mosquito_species), zi=~0, disp=~1
m_wo_temp_hum: Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature + , zi=~0, disp=~1
m_wo_temp_hum:     Relative_humidity + Host_abundance + Predator_abundance + , zi=~0, disp=~1
m_wo_temp_hum:     Blue_area_coverage + factor(Month) + (1 + Green_area_coverage | , zi=~0, disp=~1
m_wo_temp_hum:     City) + (1 | Location.Site:City) + (1 | Mosquito_species), zi=~0, disp=~1
              Df   AIC   BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
m_wo_month    15 18223 18309 -9096.4    18193                            
m_wo_temp_hum 26 18215 18364 -9081.6    18163 29.623     11   0.001815 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Removing the month worsens the fit, so our final model is:

m_final <- glmmTMB(
  Mosquito_abundance ~ Green_area_coverage * Trap_type +
                       Temperature + Relative_humidity +
                       Host_abundance + Predator_abundance +
                       factor(Month) +  
                       Blue_area_coverage +
                       (1 + Green_area_coverage | City) +          # City random intercept + slope
                      (1 | Location.Site:City) +        
                      (1 | Mosquito_species), 
  family = nbinom2(),
  data = df,
  REML = TRUE
)

summary(m_final)

 Family: nbinom2  ( log )
Formula:          Mosquito_abundance ~ Green_area_coverage * Trap_type + Temperature +  
    Relative_humidity + Host_abundance + Predator_abundance +  
    factor(Month) + Blue_area_coverage + (1 + Green_area_coverage |  
    City) + (1 | Location.Site:City) + (1 | Mosquito_species)
Data: df

      AIC       BIC    logLik -2*log(L)  df.resid 
  18332.7   18481.8   -9140.4   18280.7      2262 

Random effects:

Conditional model:
 Groups             Name                Variance  Std.Dev. Corr  
 City               (Intercept)         8.433e-02 0.290393       
                    Green_area_coverage 5.007e-05 0.007076 -0.47 
 Location.Site:City (Intercept)         9.124e-02 0.302060       
 Mosquito_species   (Intercept)         5.091e-02 0.225636       
Number of obs: 2288, groups:  City, 6; Location.Site:City, 54; Mosquito_species, 5

Dispersion parameter for nbinom2 family (): 1.98 

Conditional model:
                                         Estimate Std. Error z value Pr(>|z|)    
(Intercept)                            -0.9306526  0.3687708  -2.524 0.011614 *  
Green_area_coverage                     0.0070819  0.0036047   1.965 0.049460 *  
Trap_typeLight_CO2                      0.1241486  0.0729250   1.702 0.088677 .  
Temperature                             0.0116819  0.0082860   1.410 0.158586    
Relative_humidity                       0.0212013  0.0031738   6.680 2.39e-11 ***
Host_abundance                          0.0212067  0.0021355   9.930  < 2e-16 ***
Predator_abundance                      0.0030975  0.0017508   1.769 0.076866 .  
factor(Month)2                         -0.0609114  0.0840613  -0.725 0.468693    
factor(Month)3                          0.0934745  0.0924016   1.012 0.311724    
factor(Month)4                          0.1568507  0.1064307   1.474 0.140553    
factor(Month)5                          0.3708248  0.1202928   3.083 0.002051 ** 
factor(Month)6                          0.2913307  0.1308038   2.227 0.025932 *  
factor(Month)7                          0.4096605  0.1216106   3.369 0.000755 ***
factor(Month)8                          0.1055351  0.1059346   0.996 0.319139    
factor(Month)9                          0.1430460  0.0906455   1.578 0.114547    
factor(Month)10                         0.0811003  0.0820680   0.988 0.323051    
factor(Month)11                         0.0003518  0.0792615   0.004 0.996458    
factor(Month)12                        -0.0002625  0.0816592  -0.003 0.997435    
Blue_area_coverage                      0.0073499  0.0036179   2.032 0.042202 *  
Green_area_coverage:Trap_typeLight_CO2  0.0057406  0.0012649   4.538 5.67e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

4. Model Diagnostics

library(DHARMa)

# Simulate residuals
sim_res <- simulateResiduals(m_final)

# Plot residual diagnostics
plot(sim_res)

testDispersion(sim_res)

    DHARMa nonparametric dispersion test via sd of residuals fitted vs. simulated

data:  simulationOutput
dispersion = 0.67548, p-value = 0.728
alternative hypothesis: two.sided

testZeroInflation(sim_res)

    DHARMa zero-inflation test via comparison to expected zeros with simulation under H0 = fitted
    model

data:  simulationOutput
ratioObsSim = 1.0195, p-value = 0.912
alternative hypothesis: two.sided

ranef(m_final)

$City
            (Intercept) Green_area_coverage
Barcelona    0.05665254       -0.0006839102
Casablanca  -0.23987633       -0.0074754880
Montpellier -0.28248408        0.0087906133
Nijmegen     0.30877301       -0.0061932812
Rabat        0.06450358        0.0013067294
Utrecht      0.09243128        0.0042553368

$`Location.Site:City`
              (Intercept)
1:Barcelona    0.07043914
1:Casablanca   0.03530279
1:Montpellier -0.14612773
1:Nijmegen     0.01928858
1:Rabat        0.17883180
1:Utrecht     -0.14043348
2:Barcelona   -0.21674606
2:Casablanca   0.22542755
2:Montpellier -0.34050646
2:Nijmegen     0.52709126
2:Rabat        0.47659812
2:Utrecht      0.01212399
3:Barcelona    0.23358215
3:Casablanca  -0.17418438
3:Montpellier -0.01351448
3:Nijmegen    -0.58369434
3:Rabat       -0.06142319
3:Utrecht      0.17008047
4:Barcelona   -0.04615371
4:Casablanca  -0.01484421
4:Montpellier  0.30291202
4:Nijmegen     0.34686847
4:Rabat       -0.62301501
4:Utrecht     -0.05344528
5:Barcelona    0.09420243
5:Casablanca  -0.42233337
5:Montpellier -0.06700641
5:Nijmegen     0.18584290
5:Rabat        0.15272034
5:Utrecht      0.20916035
6:Barcelona   -0.05731574
6:Casablanca   0.16563333
6:Montpellier  0.23113679
6:Nijmegen    -0.27094925
6:Rabat        0.02313622
6:Utrecht     -0.13363323
7:Barcelona   -0.18372403
7:Casablanca  -0.48904874
7:Montpellier -0.34834010
7:Nijmegen    -0.02084802
7:Rabat       -0.31336564
7:Utrecht     -0.01039561
8:Barcelona    0.20354565
8:Casablanca   0.21546823
8:Montpellier -0.15625847
8:Nijmegen     0.22299093
8:Rabat        0.47731495
8:Utrecht      0.13005884
9:Barcelona   -0.03759367
9:Casablanca  -0.07094313
9:Montpellier  0.37974388
9:Nijmegen    -0.16354145
9:Rabat       -0.18702398
9:Utrecht      0.05690802

$Mosquito_species
                       (Intercept)
Aedes aegypti           0.38422767
Aedes albopictus       -0.14656346
Anopheles gambiae      -0.04100560
Culex pipiens          -0.16165110
Culex quinquefasciatus -0.03500751

5. Model Interpretation

library(ggeffects)
preds <- ggpredict(m_final, terms = c("Green_area_coverage", "Trap_type"))
plot(preds) + labs(y="Predicted Mosquito Abundance")

pred_temp <- ggpredict(m_final, terms = "Temperature [all]")
plot(pred_temp) + 
  labs(title = "Predicted Mosquito Abundance", 
       y="Predicted Mosquito Abundance")

pred_host <- ggpredict(m_final, terms = "Host_abundance")
plot(pred_host) +
  labs(x="Host abundance", y="Predicted mosquitoes") +
  theme_minimal()

pred_pre <- ggpredict(m_final, terms = "Predator_abundance")
plot(pred_pre) +
  labs(x="Predator abundance", y="Predicted mosquitoes") +
  theme_minimal()

pred_month <- ggpredict(m_final, terms = "Month [1:12]")

plot(pred_month) +
  scale_x_continuous(breaks = 1:12,
                     labels = month.abb) +
  labs(x="Month", y="Predicted mosquitoes",
       title="Seasonal effect") +
  theme_minimal()

NA
NA

re_plots = sjPlot::plot_model(m_final, type = "re") 

re_plots[[1]] + theme_minimal()

re_plots[[3]] + theme_minimal()

library(sjPlot)

plot_model(m_final, type = "est", 
           axis.lim = c(0.8, 1.4),        
           show.values = TRUE,             
           value.offset = 0.3,            
           value.size = 2.5,               
           dot.size = 2,                
           line.size = 0.7,
           vline.color = 'black') +              
  labs(title = "Incidence Rate Ratios (IRR)",
       x = "Incidence Rate Ratio") +
  theme_minimal()     

icc_values <- performance::icc(m_final, tolerance = 1e-6)
icc_values

# Intraclass Correlation Coefficient

    Adjusted ICC: 0.399
  Unadjusted ICC: 0.261

5. Visualizations

ggplot(df, aes(x = factor(Green_level), y = Mosquito_abundance, fill = Trap_type)) +
  geom_boxplot(alpha = 0.7, outlier.size = 0.5) +
  theme_minimal(base_size = 14) +
  labs(
    title = "Observed Mosquito Count by Greenness Level and Trap Type",
    x = "Urban Greenness Level",
    y = "Observed Mosquito Abundance"
  ) +
  scale_fill_brewer(palette = "Set2")

Mosquito abundance increases consistently with greenness level across both trap types. Although Light–CO₂ traps and Human–Odor–CO₂ traps show similar distributions, greener areas clearly host higher mosquito numbers. This provides initial, descriptive evidence of an ecological gradient before modeling.

num_df <- df %>% select(Temperature_2m, Relative_humidity, Blue_area_coverage, Green_area_coverage = Green_level)
corrplot(cor(num_df), method = "color", type = "upper", tl.col = "black")

df$Green_level <- as.factor(df$Green_level)

# Refit model
model_abundance2 <- glmmTMB(
  Mosquito_abundance ~ Green_level * Trap_type + Temperature_2m +
    Relative_humidity + Blue_area_coverage +
    (1 | Location.Site) + (1 | Mosquito_species),
  family = nbinom2(),
  data = df
)

# Estimated marginal means
emm <- emmeans(model_abundance2, ~ Green_level * Trap_type, type = "response")

# Plot interaction
plot(emm, comparison = TRUE) +
  labs(
    title = "Predicted Mosquito Abundance by Green Level and Trap Type",
    x = "Urban Greenness Level",
    y = "Predicted Abundance (Response Scale)"
  )

The model predicts a strong increase in mosquito abundance with higher urban greenness, indicating vegetation-rich environments favor mosquito populations. The effect of trap type is minimal — both traps capture similar numbers of mosquitoes across all vegetation levels. Confidence intervals widen at higher greenness levels, suggesting greater variability in mosquito counts in green zones.

# Extract random effects and standard deviations
ran_list <- ranef(model_abundance)$cond

# Convert to a tidy data frame
ran_df <- bind_rows(
  lapply(names(ran_list), function(g) {
    data.frame(
      Group = g,
      Level = rownames(ran_list[[g]]),
      Estimate = ran_list[[g]][, "(Intercept)"]
    )
  })
)

# Plot
ggplot(ran_df, aes(x = reorder(Level, Estimate), y = Estimate, color = Group)) +
  geom_point(size = 2) +
  coord_flip() +
  facet_wrap(~ Group, scales = "free_y") +
  theme_minimal(base_size = 13) +
  labs(
    title = "Random Effects by Site and Species",
    x = "Level (Site / Species)",
    y = "Deviation from Average Abundance"
  ) +
  theme(legend.position = "none",
        axis.text.y = element_text(angle = 15, vjust = 0.5),)

Weather data visualitaon

ggplot(df, aes(x = Temperature_2m)) +
  geom_histogram(binwidth = 1, fill = "steelblue", color = "white", alpha = 0.8) +
  labs(
    title = "Temperature Distribution Across Sampling Sites",
    x = "Temperature (°C)",
    y = "Frequency"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold", size = 16),
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

ggplot(df, aes(x = City, y = Relative_humidity, fill = City)) +
  geom_boxplot(alpha = 0.7) +
  labs(
    title = "Variation in Relative Humidity by City",
    x = "City",
    y = "Relative Humidity (%)"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold", size = 16),
    legend.position = "none",
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

df %>%
  ggplot(aes(x = Wind_speed_10m, fill = City)) +
  geom_density(alpha = 0.5) +
  labs(
    title = "Distribution of Wind Speed Across Cities",
    x = "Wind Speed (m/s)",
    y = "Density",
    fill = "City"
  ) +
  theme_minimal(base_size = 14)

df %>%
  mutate(Month = month(Sampling_date, label = TRUE, abbr = TRUE)) %>%
  group_by(City, Month) %>%
  summarise(Mean_Temperature = mean(Temperature_2m, na.rm = TRUE)) %>%
  ggplot(aes(x = Month, y = Mean_Temperature, color = City, group = City)) +
  geom_line(size = 1.2) +
  geom_point(size = 2) +
  labs(
    title = "Average Monthly Temperature by City",
    x = "Month",
    y = "Temperature (°C)",
    color = "City"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(face = "bold", size = 16)
  )

pred_inter <- ggpredict(
  model_abundance,
  terms = c("Green_level", "Trap_type"),
  type = "fixed",              # correct for your version
  bias_correction = TRUE
)

pred_df <- as.data.frame(pred_inter)


ggplot(pred_df, aes(
  x = factor(x),       # Greenness levels
  y = group,           # Trap types
  fill = predicted
)) +
  geom_tile(color = "white") +
  geom_text(aes(label = round(predicted, 1)), color = "black", size = 4) +
  scale_fill_gradient(low = "#fefefe", high = "#01665e") +
  labs(
    title = "Predicted Mosquito Count by Greenness Level and Trap Type",
    x = "Urban Greenness Level",
    y = "Trap Type",
    fill = "Predicted\nAbundance"
  ) +
  theme_minimal(base_size = 14)+
  theme(plot.title = element_text(size=12))

Mosquito abundance rises markedly with increasing urban greenness, regardless of the trap type used. The model predicts roughly five times more mosquitoes in highly vegetated areas than in low-greenness zones, highlighting the ecological importance of vegetation in mosquito distribution.

df_pred <- df %>%
  mutate(
    Predicted = predict(model_abundance, type = "response"),
    Observed = Mosquito_abundance
  )


ggplot(df_pred, aes(x = Predicted, y = Observed)) +
  geom_point(alpha = 0.4, size = 2) +
  geom_smooth(method = "lm", color = "#01665e", se = TRUE) +
  labs(
    title = "Model Fit: Predicted vs Observed Mosquito Abundance",
    x = "Predicted Abundance (Model)",
    y = "Observed Abundance (Data)"
  ) +
  theme_minimal(base_size = 14)

df %>%
  group_by(Green_level, Mosquito_species) %>%
  summarise(Mean_abundance = mean(Mosquito_abundance)) %>%
  ggplot(aes(x = factor(Green_level), y = Mean_abundance, fill = Mosquito_species)) +
  geom_bar(stat = "identity", position = "fill") +
  labs(
    title = "Relative Mosquito Species Composition by Greenness Level",
    x = "Urban Greenness Level",
    y = "Proportion of Total Abundance"
  ) +
  theme_minimal(base_size = 14)

df %>%
  mutate(Month = month(Sampling_date, label = TRUE, abbr = TRUE)) %>%
  group_by(Month) %>%
  summarise(Mean_abundance = mean(Mosquito_abundance)) %>%
  ggplot(aes(x = Month, y = Mean_abundance, group = 1)) +
  geom_line(size = 1.2, color = "#01665e") +
  geom_point(size = 2, color = "#01665e") +
  labs(
    title = "Seasonal Trend in Mosquito Abundance",
    x = "Month",
    y = "Average Abundance"
  ) +
  theme_minimal(base_size = 14)

df %>%
  mutate(Month = month(Sampling_date, label = TRUE, abbr = TRUE)) %>%
  group_by(City, Month) %>%
  summarise(Mean_abundance = mean(Mosquito_abundance, na.rm = TRUE)) %>%
  ggplot(aes(x = Month, y = Mean_abundance, group = 1)) +
  geom_line(color = "#01665e", size = 1.2) +
  geom_point(color = "#01665e", size = 2) +
  facet_wrap(~ City, ncol = 2) +
  labs(
    title = "Seasonal Trends in Mosquito Abundance by City",
    x = "Month",
    y = "Average Mosquito Abundance"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold", size = 16),
    strip.text = element_text(size = 12),
    axis.text.x = element_text(angle = 90, vjust = 0.5)
  )

df %>%
  mutate(Month = month(Sampling_date, label = TRUE, abbr = TRUE)) %>%
  group_by(City, Month) %>%
  summarise(Mean_abundance = mean(Mosquito_abundance, na.rm = TRUE)) %>%
  ggplot(aes(x = Month, y = Mean_abundance, color = City, group = City)) +
  geom_line(size = 1.2) +
  geom_point(size = 2) +
  labs(
    title = "Seasonal Trends in Mosquito Abundance by City",
    x = "Month",
    y = "Average Mosquito Abundance",
    color = "City"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold", size = 16),
    legend.position = "bottom"
  )

df %>%
  mutate(Month = month(Sampling_date, label = TRUE, abbr = TRUE)) %>%
  group_by(Mosquito_species, Month) %>%
  summarise(Mean_abundance = mean(Mosquito_abundance, na.rm = TRUE), .groups = "drop") %>%
  ggplot(aes(x = Month, y = Mean_abundance, color = Mosquito_species, group = Mosquito_species)) +
  geom_point(size = 2) +
  geom_smooth(se = FALSE, linewidth = 1.2, method = "loess") +
  labs(
    title = "Seasonal Trends in Mosquito Abundance by Species",
    x = "Month",
    y = "Average Mosquito Abundance",
    color = "Species"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold", size = 16),
    legend.position = "bottom"
  )

df$Date <- as.Date(df$Sampling_date, format = "%Y-%m-%d")

# Create the monthly trend plot
df %>%
  mutate(Month = month(Date, label = TRUE, abbr = TRUE)) %>%
  group_by(Mosquito_species, Month) %>%
  summarise(Mean_abundance = mean(Mosquito_abundance, na.rm = TRUE), .groups = "drop") %>%
  ggplot(aes(x = Month, y = Mean_abundance, color = Mosquito_species, group = Mosquito_species)) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 2) +
  labs(
    title = "Monthly Trends in Mosquito Abundance by Species",
    x = "Month",
    y = "Average Mosquito Abundance",
    color = "Species"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold", size = 16),
    legend.position = "bottom"
  )

# Create faceted plot by species
df %>%
  mutate(Month = month(Date, label = TRUE, abbr = TRUE)) %>%
  group_by(Mosquito_species, Month) %>%
  summarise(Mean_abundance = mean(Mosquito_abundance, na.rm = TRUE), .groups = "drop") %>%
  ggplot(aes(x = Month, y = Mean_abundance, group = 1)) +
  geom_line(color = "#2C7BB6", linewidth = 1.2) +
  geom_point(color = "#2C7BB6", size = 2) +
  facet_wrap(~ Mosquito_species, scales = "free_y") +
  labs(
    title = "Monthly Trends in Mosquito Abundance by Species",
    x = "Month",
    y = "Average Mosquito Abundance"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
    strip.text = element_text(face = "bold"),
    axis.text.x = element_text(angle = 90, vjust = 0.5),
    legend.position = "none"
  )

ggplot(df, aes(x = Trap_type, y = Mosquito_abundance, fill = Trap_type)) +
  geom_boxplot(outlier.alpha = 0.2) +
  facet_wrap(~ City) +
  scale_fill_brewer(palette = "Set2") +
  labs(
    title = "Distribution of Mosquito Abundance by Trap Type and City",
    x = "Trap Type", y = "Mosquito Abundance"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    legend.position = "none",
    axis.text.x = element_text(angle = 45, hjust = 1)
  )